A new constant-pressure molecular dynamics method for finite system 
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In this letter, by writing the volume as a function of coordinates of atoms, we present a new constant- 
pressure molecular dynamics method with parameters free. This method is specially appropriate for 
the finite system in which the periodic boundary condition does not exist. Simulations on the 
carbon nanotube and the Ni nanoparticle clearly demonstrate the validity of the method. By using 
this method, one can easily obtain the equation of states for the finite system under the external 
pressure. 
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The molecular dynamics (MD) simulation method is 
a powerful tool, widely used in chemistry, physics, and 
materials science. As a very important achievement, 
constant-pressure MD proposed by Andersen, and sub- 
sequently extended by Parrinello and Rahman, ^ has 
provided a powerful method to study systems under the 
external pressure. Now the constant-pressure MD has 
become a standard tool to study the physical properties 
under the external pressure. Especially, it has played a 
central role in investigating the structural phase transi- 
tion as well as in geophysical and astrophysical applica- 
tions. 

Recently, the study of low-dimensional and biological 
system under the external pressure become a considerable 
interest. |^,|^ In particular, studies on the carbon nan- 
otube, H] clusters and nanocrystal, such as CdSe, CdS 
and Si nanocrsystals, 0] have revealed a wealth of in- 
teresting new phenomena. Usually computer simulation 
could substantially complement the experimental infor- 
mation. However, the traditional constant-pressure MD 
method was designed for the infinite system with a pe- 
riodic boundary condition, it can not be directly applied 
to the finite system, such as nanotubes and nanocrystals, 
in which the periodic boundary condition does not exist. 
For this reason it is necessary to develop a new computa- 
tional scheme to study the finite system. Very recently, 
Martonak, Molteni and Parrinello have successfully stud- 
ied the pressure induced amorphorpization of Sis^Hse 
cluster by introducing pressure-transmitting liquid. ||] In 
order to have a well-defined isotropic constant pressure on 
the cluster, the number of particles and the volume of the 
pressure-transmitting liquid should be much larger than 
that of the cluster, thus it costs a significant overhead of 
the computation. Additionally, one also must determine 
how the liquid atoms interacts with themselves and with 
cluster atoms so that the liquid does not crystallize, or, 
in the time scale of the simulation the liquid does not 



undergo a glass transition. 

In this paper, by writing the volume as a function of 
coordinates of atoms, we propose a new constant-pressure 
MD, which should be appropriate for finite systems with 
arbitrary shape. This new constant-pressure MD is pa- 
rameter free, and can be used in any system with arbi- 
trary shape, especially in nanocrystals. We will demon- 
strate its validity by classical simulations of a carbon nan- 
otube, and small Ni particles, its extension to ab — initio 
molecular dynamics method is straight forward. Q 

We write the Lagrangian L of a TV-atom system as: 
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L = y^ 
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where rj, m,, Pi are the coordinates, mass, momentum 
of ith atom respectively, and is potential of system. 
Y is the volume of the system, and P^xt is the external 
pressure. 

If the system obeys the Newton's mechanism, the equa- 
tions of motion for derived from the Lagrangian L read. 



(2) 



Obviously the enthalpy will be conserved. The equa- 
tions of motion derived from Eq. 2 produce the constant 
pressure ensemble for the system, as we show bellow. For 
an equilibrium system with the external pressure, 

^ N N N 

— m,v^ - ^ r, • V0 - E • ^-tV^) = (3) 
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Then we have, 
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It is well known that, in statistical physics, the 
volume is an additive quantities, which can be written as 
a summation of the volume of individual atoms, 

N 
i 

where Vi is the volume of the ith atom, which can be 
generally written as a cubic homogeneous function of its 
nearest neighbor distance r^ , i.e. 

According to Euler theorem, 

N 

J2r,-vV^3V (7) 
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So from Eq. 4 we have 

^ N N 

Pext = Pint = ^ (E '"'^'^ -H^^- V<?^) (8) 
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where Pint refers to the internal pressure, since the ex- 
ternal pressure Pext is a constant, Pint is also a constant. 
Thus, by writing the volume as a function of atomic co- 
ordinates, a new constant-pressure MD is presented, in 
which no extra parameter is introduced. 

The present constant-pressure MD method has several 
advantages. First, this method can make the calculation 
more realistic, without needing to choose the mass for the 
volume as in traditional constant-pressure MD, |^ which 
directly affects the time scale of the relaxation. Secondly, 
in the present method, the responding of system to the 
external pressure is more physical. This is specially im- 
portant for the inhomgenous system. However, in the 
traditional constant-pressure MD, the responding of sys- 
tem to the external pressure is essentially linear, i.e., the 
volume of the system is linearly scaled according to the 
difference of the internal and external pressure. 

The key to the success of the present constant-pressure 
MD method is to define the volume as a function of 
atomic coordinates. In the traditional constant-pressure 
MD method, the volume is a generalized coordinates, 
which has equal importance as an atomic coordinate, 
the constant pressure is dynamically achieved by directly 
changing the volume of system. However, in the present 
scheme, the constant pressure is obtained by dynamically 
changing the motion of each atom. 

To show how well this new constant-pressure MD works 
in the real application, we have simulated carbon nan- 
otubes and Ni nanocrystals under the external pressure. 

In the real simulation, to get an exact formalism for the 
atomic volume is physically difficult. However, there are 
a varity of sufficient approximations. One of the simplest 



and direct way is based on the Wigner-Seitz premitive 
cell, we used the scaled volume of the atomic sphere to 
approximate the Wigner-Seitz premitive cell, which has 
the following form, 

where the summation runs over all the first nearest neigh- 
bors of the zth atom, ji is a scale factor which is closely 
related the number of the nearest neighbours of the zth 
atom. At the high temperature or in the system where 
the first nearest neighbors frequently change, we used a 
more sophisticated form. |^ 

The definition of the volume is not unique, our defini- 
tion for Ni nanocrystals is one of them. If it is necessary, 
one can make more elegant one. In the study of carbon 
nanotube, we fix the periodic length in axial direction. 
The slight different definition of the volume for nanotubes 
will be presented else. Q 

The interaction between carbon atoms is described by 
the parametrized potential developed by Brenner ac- 
cording to the Tersoff bonding formalism which is 
widely used to study mechanical properties of carbon 
nanotube [Q. The Sutton-Chen potential is used to de- 
scribe the interaction between Ni atoms, which is also 
widely used in the literature jl^. We find that, the ob- 
tained results for the carbon nanotube and Ni nanocrys- 
tal are in good agreement with what obtained by other 
methods. 

The evolution of the instantaneous pressure in the 
simulation of the (10x10) carbon nanotube and the Ni 
nanocrystal (Fig. 1), shows that the present method 
does recover a constant pressure simulation. Although 
the instantaneous pressure of the system fluctuates, their 
average values are equal to the set external pressure, 0.7 
GPa for the carbon nanotube and 7.0 GPa for the Ni 
nanocrystal. The correlation between the volume and 
pressure can also be clearly observed. 

What shown in Fig. 2 is the calculated energy via vol- 
ume (pressure) for carbon nanotube, the close agreement 
between the static calculation and the present constant- 
pressure MD simulation clearly demonstrates the validity 
of the present method. In the static calculation, energies 
are calculated with linear scaling the radius of the carbon 
nanotube without relaxing the atomic positions. All the 
energies shown in the lower panel of Fig. 2 are relative to 
the minimum energies, and the volume are renormalized 
by the equilibrium volume without the external pressure. 
We performed the constant-pressure MD simulation at 
the zero temperature and 300 K at various external pres- 
sure from GPa to 2 GPa. Both the static calculation 
and the constant-pressure MD simulation at the zero tem- 
perature give the very similar results, the slight difference 
comes from the relaxation of the atomic coordinates. 

We show the enthalpy as a function of the reduced 



volume for static and present MD results in the middle 
panel of Fig. 2. It can be seen that, at each volume, the 
enthalpy of present MD, in which the structures under 
the external pressure are relaxed, is always smaller than 
that from the static calculation as it should be. The 
difference of the calculated enthalpy between the static 
and MD at K increasing with reducing volume suggests 
the atomic relaxation becomes more and more important. 
The reasonable agreement between the ab-initio and the 
model potential calculation shows that the present model 
potential can qualitatively describe the behavior of the 
carbon nanotube. The new constant-pressure MD also 
correctly describe the finite temperature properties. 

The equation of states (EOS) at finite temperature for 
Ni nanocrystal can be obtained through the simulation, 
which is reasonable in agreement with the bulk phase 
EOS calculated by traditional constant-pressure MD. Fig. 
3 shows the equation of states for a Ni nanocrystal and 
bulk phase calculated by the new constant-pressure MD 
and traditional MD respectively. All the energies shown 
in Fig. 3 are relative to the minimum energies, and the 
volume is renormalized by the equilibrium volume with- 
out the external pressure at 300 K. The nanocrystal and 
bulk phase show similar behavior. However, from the fig- 
ure, one can see that the nanocrystal is not as hard as 
the bulk phase. If assuming first-order Birch-Murnaghan 
EOS. ijl^ is also appUicable to the finite system, we yields 
the bulk modulus 136 GPa and 166 GPa for nanocrystal 
and bulk phase respectively, in agreement with our pre- 
vious results. ||l^ 

In summary, by writing the volume as a function of 
coordinates of atoms, we have proposed a new constant- 
pressure molecular dynamics method for the finite system 
with parameters free, where the external pressure could 
be exactly implemented. Simulations on the carbon nan- 
otube and Ni nanoparticles have clearly demonstrated the 
validity of the method, in which the constant pressure is 
recovered. We have also shown that, with the newly pro- 
posed scheme, the EOS for the finite system can be cal- 
culated through a molecular dynamics simulation. The 
extension to ab — initio molecular dynamics method is 
straight forward and is in progress. 

We thanks Dr. D. J. Shu and Dr. G. Chen for the 
technical assistance. This work is supported by the Na- 
tional Science Foundation of China, the special funds for 
major state basic research and CAS projects. 
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FIG. 1. Tlie evolution of the instantaneous pressure and 
volume through MD runs for the carbon nanotube (lower two) 
and Ni nanocrystal (upper two). The pressure and volume 
of the system fluctuates around the average value. The new 
constant-pressure MD does recover a constant pressure. 
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FIG. 2. The calculated properties for the carbon nanotube. 
Low panel: The energy as a function of the volume by the 
new constant-pressure MD at K (filled circle), 300 K (open 
square), static calculation (open circle) and ab — initio re- 
sult (solid line). Middle panel: The enthalpy as a function of 
the reduced volume for static (open circle) and present MD 
(filled circle) calculation. Up panel: The pressure-volume re- 
lation of the carbon nanotube at K (filled circle) and 300 
K (open square). Our MD results at zero temperature are in 
good agreement with the static results except for large vol- 
ume charge, where the structural relaxation includes in our 
MD runs. 
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FIG. 3. The equation of states for Ni nanocrystals (filled 

circles) and bulk phase (open circles), where the data for 
Ni-jir,! nanocrystals and bulk phase are calculated by the new 
and traditional molecular dynamics simulation respectively. 
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